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We report the results of a molecular dynamics study of the ionic liquid l-#-butyl-3-m^ 
hexafluorophosphate [>mim][PF 6 ], a widely studied ionic liquid. An all-atom force field is developed using 
a combination of density functional theory calculations and CHAKMM 22 parameter values. Molecular 
dynamics simulations are carried out in the isothermal— isobaric ensemble at three different temperatures. 
Quantities computed include infrared frequencies, molar volumes, volume expansivities, isothermal com- 
pressibililties, self-diffusivities, cation— anion exchange rates, rotational dynamics, and radial distribution 
functions. Computed thermodynamic properties are in good agreement with available experimental values. 



1. Introduction 



The term l ionic liquidl (IL) has been coined in recent years to 
describe a class of organic salts that are liquid in .their pure 
state at or near room temperature. 1 Some of the more widely 
studied ILs consist of a heterocyclic cation based on a substituted 
pyridine or imidazole and an inorganic anion. Early studies 
focused on compounds with the [AlCU]" anion, but these liquids 
proved to be unstable in air. This problem was- overcome 2 
through the use of alternative anions such as [BF 4 ]~ [N0 3 ]~, 
[CH3COO]", and [PF 6 r. These stable ILs have attracted a great 
deal of interest because of their potential as nonvolatile (and 
hence potentially environmentally benign) solvents. 3 They also 
show potential for a range of other applications including 
separations, industrial cleaning, fuel cells, lubricants, and heat 
transfer fluids. 4 

Compared to conventional organic solvents, relatively little 
is known about the thermodynamic and transport properties of 
ILs and how these properties relate to chemical constitution and 
structure. In addition to experimental efforts directed- at this 
problem, a few theoretical studies have been reported in which 
quantum chemical 5 and classical condensed phase simulations 6 " 8 
were used to examine structural and thermophysical properties 
of water-stable ILs. Other theoretical studies have focused on 
chloroaluminate-based ILs. 9-11 

The present work reports results of single molecule quantum 
chemical calculations and condensed phase classical molecular 
dynamics (MD) simulations on the ionic liquid l-n-butyl-3- 
memylimidazolium hexafluorophosphate, which will be abbre- 
viated as [bmim][PF6]. Figure 1 shows the optimized geometry 
of a single [bmimlfPFs] molecule in the gas phase, obtained 
from a density functional theory calculation (see below). The 
atom labels shown in this figure will be referred to throughout 
the work. 

2. Force Field 

■The extent to which a classical molecular simulation ac- 
curately predicts thermophysical properties depends on the 
quality of the force field used to model the interactions in the 
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Figure 1. Geometry optimized structure of J>mim][PF6] with atom 
labels noted. 

fluid. In this work, we have used a standard molecular 
mechanics force field, 12 with functional form 
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where V tot is the total energy of the system, harmonic potentials 
govern bond length, bond angle, and improper angle motion 
about nominal values r 0 , 6q, and ipo and dihedral angles are 
modeled using, a standard cosine series. The Lennard- Jones 
parameters for unlike atoms, 6,y and r m i n .y, are obtained using 
the Lorentz— Berthelot combining rule. Coulombic interactions 
are modeled using fixed partial charges on each atom center. 



10.1021/jp0267003 CCC: $22.00 © 2002 American Chemical Society 
Published on Web 11/16/2002 



BEST AVAILABLE COPY 



12808, /. Phys. Chem. B, Vol. 106, No. 49, 2002 



Morrow and Maginn 



TABLE 1: Partial Atomic Charges and Lennard- Jones 
Parameters Used in This Work 
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With the exception of ro, #o> and #0 all cation intramolecular 
and Lennard- Jones parameters in eq 1 were taken directly from 
the CHARMM 22 force field 13 and used without modification. 
For the anion, bond length and bond angle parameters were 
derived from ab initio calculations, as described below. 
CHARMM 22 Lennard- Jones parameters were used for phos- J 
phorus and fluorine. 

The minimum-energy geometry of the |>mim] cation and the 
\?F$] anion was determined by performing ab initio geometry 

TABLE 2: Bond, Angle, Dihedral, and Improper Force Constants 



optimizations on the isolated cation and anion at the B3LYP/ 
6-31 1+G* level of>theory using Gaussian 98. 14 Cation and anion 
charges were set at +1 and -1, respectively. The resulting 
imnimized structure was used to set ro, 6o, and ip 0 . Geometry 
optimization was also carried out on the [bmimJCPFe] pair at 
the same level of theory. The anion was initially located weD 
away (over 10 A) from the cation and facing the imidazolium 
C 2 -carbon side. The rninimum energy structure agrees very well 
with the results of Meng et a!. 5 Partial atomic charges were 
derived from this geometry using the CHELPG 15 method. A 
listing of all force field parameters is given in Tables 1 and 2. 
We note that the pair calculations yielded total charges on the 
cation and anion of +0.904 and —0.904, respectively. In 
addition, the anion is slightly polarized (i.e., the charges on the 
fluorine atoms are not symmetric). For most of the simulations, 
we used these partial charges as a simple approximation of the 
induced polarization of the ions that occur in the liquid. This 
turned out not to be essential, given that the anion exhibited no 
preferential orientation relative to the nearest cation, as deter- 
mined by monitoring the fraction of the time a given fluorine 
atom was closest to the C2 carbon of the cation. We also 
performed simulations on an anion with symmetric charges 
totaling -0.904 and found essentially no difference in the 
computed properties. 

The force constants kt, and k& for the anion were determined 
by performing a vibrational analysis on the geometry-optimized 
pair within Gaussian 98. The normal modes corresponding to 
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the stretching of a P-F bond and the bending of an F-P— F 
bond were identified using Gaussview, a graphical interface for 
the Gaussian program, and the computed vibrational frequencies 
were used to calculate the corresponding harmonic force 
constants. The geometry and force constants for the cation and 
anion are summarized in Table 2. In a similar manner, 
vibrational assignments for the major fundamental modes of 
the pair were made and compared against experimental infrared 
spectroscopy measurements, 16 - 17 as described in the Results . 
section. 

£ 3. Simulation Details 

Molecular dynamics simulations of 300 [brnim] cations and 
300 [PF6] anions were performed with the program NAMD 18 
version 2.4 in a cubic cell with standard periodic boundary 
conditions. The simulations were carried out in the isothermal— 
isobaric (NPT) ensemble using a Nose-Hoover barostat to 
maintain a pressure of 0.98 atm. The temperature was controlled 
via Langevin dynamics with a damping factor of 5 ps -1 . Initial 
configurations were generated by randomly inserting cation— 
anion pairs into the simulation box. To prevent overlap, an 
insertion was rejected if any of the newly inserted atoms were 
within 0.75 A of any existing atoms. The initial configurations 
were relaxed using a conjugant— gradient energy minimization 
scheme. For most of the simulations, the initial cell volume was 
chosen to match the experimental density of the state point of 
interest. All of the C— H bonds were held rigid using the 
SHAKE 19 algorithm. The r-RESPA 20 multiple time-stepping 
algorithm was used with a time step of 2 fs for bonded and van 
der Waals interactions and 4 fs for electrostatic interactions. 
The dispersion interactions were cut off beyond 20.0 A. A 
switching function, initiated at a distance of 18.5 A, was used 
to bring the dispersion interactions smoothly to zero at the cutoff 
distance. Long-range electrostatic interactions were computed 
using the particle mesh Ewald method. 21,22 

Prior to conducting a production run, the total energy and 
cell volumes were monitored until steady-state was reached. 
Equilibration times varied from 700 to 1000 ps, after which 
production runs lasting 4 ns were started. The thermodynamic 
properties of the system (total energy, pressure, temperature, 
kinetic, and potential energy contributions) were saved to disk 
every 100 time steps, and the atomic coordinates were saved to 
disk every 400 time steps. All classical simulations of [bmim]- 
[PF6] were conducted at atmospheric pressure (0.98 bar) and at 
the temperatures 298, 323, and 343 K. 

The program NAMD was chosen primarily because it is a 
parallel MD program that scales remarkably well with the 
number of processors. Simulations were run on a cluster of eight 
Dell Precision 530 computers running Linux. Each node contains 
two Intel Xeon processors operating at 1.7 GHz. A 4 ns 
simulation of 300 IL molecules (9600 atoms) took approxi- 
mately 300 h (12.5 days) to complete. 

4. Results and Discussion 

Figure 2 shows the experimental 16,17 and computed IR spectra 
for [bmim] [PFg]. The relevant peaks in the spectra are labeled 
I— V. The computed frequencies of the major vibrational modes 
agree well with experiment. Using Gaussview, peak I was 
visually identified as the out-of-plane bending of the C2— Nj — 
C5 angle in the imidazolium ring, peak II is the F-P— F bending 
motion, peak III is the in-plane bending of H-C— N angles in 
the imidazolium ring, peak IV is the vibrational motion of 
several atoms in the imidazolium ring, and the peaks in region 
V are the "stretching of C— H bonds. The computed frequency 
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Figure 2. Comparison of computed and experimental 16 IR. spectra for 
[bmim][PF 6 ]. 

TABLE 3: Comparison of Predicted and Experimental IR 
Spectral Data 

pred. freq. expt. 16 freq. 
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Figure 3. Comparison of predicted and experimental molar volumes 
of (>rnim][PF6] at 0.98 bar. Lines are linear fits to the data. 

for peak I occurs at 50 cm" 1 less than the experimental value. 
The agreement between computed and experimental frequencies 
for peaks n— IV is excellent. Peaks II and III occur at 15 cm" 1 
less. than the experimental values, and peak IV occurs at 26 
cm" 1 less than the experimental value. The computed frequen- 
cies above 2800 cm" 1 (region V) are overestimated, with errors 
larger than 100 cm" 1 (see Table 3). The large errors in the 
predicted C— H stretching frequencies are likely due to anhar- 
monicity effects. 23 

4.1. Molar Volume. Figure 3 shows the experimental 24 and 
computed molar volumes as a function of temperature at 
atmospheric pressure. Results are also given in Table 4. At all 
temperatures, the predicted molar volumes are lower than the 
experimental molar volumes by less than 1%. This level of 
agreement is excellent considering that the calculations are 
purely predictive; none of the force field' parameters have been 
adjusted to match experimental data. Table 5 gives the 
breakdown of the potential energy into electrostatic, van der 
Waals, intramolecular, and kinetic energies. The average 



12810 J. Phys. Chem. B, Vol 106, No. 49, 2002 

TABLE 4: Comparison of Predicted and .Experimental 
Molar Volumes. Volume Expansivities, and Isothermal 
Compressibilities at 0.98 Bar 
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TABLE 5: Comparison of the Contribution of Various 
Terms to the System Total Energy 



temp 

(K) 



elect. 
(kJ/mol) 



LJ 
(kJ/mol) 



kinetic 
(kJ/mol) 



intra. 
(kJ/mol) 



298.2 
323.2 
343.2 



-256.72 
-256.21 
-255.07 



-89.22 
-87.19 
-85.51 



100.50 
108.92 
116.78 
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102.4 



electrostatic contribution to the potential energy at 298 K is 
—256.7 kJ/mol, which is twice as large as the van der Waals 
interactions and is 46 kJ/mol larger than the value reported for 
[ernimjfAlCU]. 9 The total intramolecular potential energy, 
calculated by surrrrning the bond, angle, dihedral, and improper 
angle energies, is 89.6 kJ/mol. This suggests that properties are 
most sensitive to the electrostatic portion of the force field. 

4.2. Volume Expansivity (ctp). The volume expansivity 
quantifies the extent to which the volume of a fluid changes 
with temperature at constant pressure, and is defined as 



1(8V\ 



(2) 



The volume expansivity can. be computed by running a series 
of simulations at the same pressure but different temperatures-. 
Because the change of volume with temperature is approxi- 
mately linear, ctp can be calculated using eq 2 by fitting a straight 
line to the simulated molar volume data. Table 4 lists the volume 
expansivities computed using this approach compared with the 
experimental results. 24 In all cases, the simulations predict 
expansivities that are lower than the experimental values by 
about 0.6 K" 1 . It is observed that the predicted expansivities 
decrease slightly with increasing temperature, in good agreement 
with the experiments. 

4.3. Isothermal Compressibility (/Cx). The isothermal com- 
pressibility quantifies the extent to which the volume of a fluid 
changes with pressure at constant temperature, and is defined 
as 



(3) 



In this work, the isothermal compressibility was calculated using 
the following fluctuation formula 25 



(4) 



Table 4 lists the computed isothermal compressibilities com- 
pared with experimental values. 24 The computed isothermal 
compressibilities are consistently lower than the experimental 
values by 12-33%. There is still reasonably good agreement 
with experiment, however, especially considering the well- 
known difficulty in computing pressures from an atomistic 
simulation, as well as the inherent inaccuracy involved in 
computing derivative quantities using a fluctuation formula of 
the type in eq 4. The overall good agreement between the 
calculated and experimental PVT properties for this fluid 
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TABLE 6: Self-Diffusion Coefficients and Cohesive Energy 
Density of PbmimjfPFe] versus Temperature 
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Figure 4. Center of mass radial distribution functions of [bmimJCPFe] 
at 298 K and 0.98 bar. 

indicates that the force field does a reasonable job of describing 
the liquid state of [brnimJfPFe]. The force field was therefore 
used to compute other thermodynamic, structural, and dynamic 
properties for which experimental data do not yet exist. 

4.4. Cohesive Energy Density. The cohesive energy density 
of a liquid is defined as 26 



c — ■ 



A£/ vap 



(5) 



where AlT* p is the change in internal energy of the liquid upon 
isothermal vaporization into the ideal gas state and V 1 - is the 
molar volume of the liquid. To calculate the cohesive energy 
density, the average internal energy and molar volume of the 
liquid were obtained from the liquid simulation, and the average 
internal energy of the ideal gas was obtained by simulating a 
single [bmim][PF6] ion pair at the same temperature as the liquid 
but at zero pressure (i.e., no periodic boundary conditions). The 
cohesive energy density of [bmim][PF6] versus temperature is 
given in Table 6. At 298 K, the cohesive energy density is 761 
J cm -3 . In contrast, the cohesive energy densities of the heavy 
hydrocarbons hexadecane and naphthalene are 268- and 410 J 
cm" 3 , respectively. 27 The extremely high cohesive energy 
density of [^mimJtPFe] stems mainly from electrostatic interac- 
tions, and explains why these liquids have such low volatility. 

4.5. Molecular Structure. To obtain a better understanding 
of the structure of this ionic liquid, various radial distribution 
functions (RDFs or g(r)) were computed at each of the 
statepoints. The center-of-mass RDFs for cation— cation, cation— 
anion, and anion— anion pairs at 298 K and 0.98 bar are shown 
in Figure 4. The RDFs at 323 and 343 K are qualitatively very 
similar and are not shown here. It is observed that the first 
solvation shell for cation— anion pairs forms at about 4.3 A. 
The second and third cation— anion solvation shells occur at 
10.6 and 17^6 A, respectively. Given the long-range Couiombic 
interactions in this system, it is - not surprising to see that weak 
ordering persists beyond 23 A. The first peak in the cation- 
cation RDF occurs at 8 A, and the anion— anion RDF shows 
two peaks at 6.5 and 9 A. This split peak is the result of 
sequential ordering induced by the cation-anion pairs. The 
RDFs agree remarkably well with those of Shah et al., 8 who 
used a simpler, united atom model for this liquid. 



Dynamics Study of [bmimjfPFg] 



/. Phys. Chem. B, Vol 106, No. 49, 2002 12811 



-s 3 

2 
1 
0 





'"!' 










P-C 4 

....... P-C 5 

P-C a 








L_ 














































































• i 
r 



















^.0 2 4 6 8 10 12 14 16 18 20' 
r(A) 

Figure 5. Atom-atom radial distribution functions for [brnim][PF 6 ] 
at 298 K and 0.98 bar. 

Further insight into the liquid structure can be gained by 
examining the site- site pair RDFs. RDFs for the phosphorus 
atom of the anion and different imidazolium carbon atoms of 
the cation are shown in Figure 5. It is observed that the anion 
prefers to associate with the C2 carbon of the imidazolium ring. 
The C2-H1 atom pair constitutes the largest amount of positive 
charge on the cation, so it is not surprising that the anion 
interacts more strongly with this part of the cation. In addition, 
the locations of the first peaks in the RDFs of Figure 5 are 
consistent with the ab initio optimized geometry of the cation - 
anion pair shown in Figure 1 . 

The RDFs suggest that the individual ions in [bmim][PF6] 
tend to aggregate into well-defined ion clusters. By integrating 
g(r) out to the location of the first minimum, the coordination 
number for the first solvation shell, N, can be calculated via 
the following equation 



N = f Q rtiKU pg(r)4jzr 2 dr 



(6) 



where p is the bulk density and r S h e ii is the first minimum in 
g(r). Using 7- she u equal to 8.36 A, the calculated coordination 
number for the cation— anion first solvation shell is 6.8, 
indicating that each ion is surrounded by a cage of nearly seven 
■ other counterions. Similarly, the cation— cation coordination 
number is 6.2, and the anion— anion coordination number is 6.1. 
Therefore, the total coordination number for the first solvation 
shell of an ion is nearly 13. Note, however, that the coordination 
numbers computed from eq 6 are sensitive to the choice of r she ii- 
For example, setting r shc ii equal to 8.0 A reduces the cation— 
anion, cation ^cation, and anion— anion coordination numbers 
to 6.2, 5.0, and 4.9, respectively. 

4.6. Diffusion. The coefficient of self- diffusion for a fluid 
or solid can be calculated using the Einstein relation 28 



D Setf = |lim|(|r ( -(0-^(0)l> 



(7) 



where the quantity in braces is the ensemble- averaged mean 
square displacement (MSD) of the molecules and r, is the vector 
coordinate of the center of mass of ion i. The mean square 
displacements of the cation and anion at 298 K and 0.98 bar 
are shown in Figure 6. The MSDs are linear up to" roughly 1 
ns, after which they. show a slightly nonlinear increase- up to 
about 1.8 ns, followed by a second linear regime. A block 
averaging technique 25 was used to obtain the MSDs,. so the 
results are most reliable at short times. For this reason, the self- 
diffusion coefficients for the cation and anion were obtained 
by fitting the slope of the linear region from 200 to 1000 ps to 
a straight line and applying eq 7. Results are given in Table 6. 
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Figure 6. Center of mass mean square displacements of [brnim] [PF 6 ] 
at 298 K and 0.98 bar. 

At 298 K and 0.98 bar, the predicted self-diffusion coefficients 
of the cation and anion are 9.70 x 10" 12 and 8.82 x 10" 12 m 2 
s" 1 , respectively. For comparison, the self-diffusion coefficient 
of water 29 at 298 K and 1 atm is 2.3 x 10" 9 m 2 s** 1 . This result 
is consistent with the fact that [bmim][PF6] has a much higher 
viscosity (450 cP) 30 than water. If one assumes that the Stokes— 
Einstein relation can be applied to this system, then the self- 
diffusion coefficients for the ionic liquid can be estimated from 
the relation. 



(8) 



where Djl = l/2(Z) [b miin] + Apf 6 j) and V * s ^ e viscosity. Using 
eq 8 and rju 2 o = 0.9 cP, 31 the estimated self-diffusion coefficient 
for the ionic liquid is 4.6 x 10" 12 m 2 s~ J , which is fairly close 
to the computed value, especially considering the simphfying 
assumptions behind the Stokes— Einstein relation. The computed 
self-diffusion coefficients for [brnim] and [PF 6 ] are roughly 10 
times smaller than those reported by de Andrade et al., 9 who 
computed self-diffusivities for l-emyl-3-memylimidazolim 
([emim]) and [AlCU] ions and are also 10 times smaller than 
the experimental result for [emim]. 32 This is not surprising 
because the [emim] and [AICI4] ions are smaller than the ions 
in this. study and because the reported intermolecular potential 
energy 9 for [eminrHAlCLJ is not as great as it is in this system. 
The computed self-diffusion coefficient for [PF6] is also seven 
times lower than that calculated by Hanke et al. 6 from a 
simulation of dimemylimidazolium ([dmim]) [PFe] at 400 K. 
Much of this difference could be due to the fact that the 
simulations of Hanke et al. 6 were at higher temperatures than 
the present work and because [dmim] is a smaller cation than 
[brnim]. We also note that the simulations in both ref 6 and 9 
were over a much shorter time scale (100 ps) compared to the 
- relatively long simulation times of this work. 
■ The self-diffusion coefficients for [>iriim] and [PF6] could 
not be determined reliably by fitting the mean square displace- 
ment data at times longer than 1 ns. This was due to inaccuracies 
in the data caused by insufficient sampling. These inaccuracies 
were apparent in the relatively large (but not systematic) 
anisotropy in the diagonal terms of the self-diffusivity tensor 
at times greater than 1 ns. Long-time anisotropy in a homoge- 
neous system indicates that even longer simulations than those 
in this work are needed to accurately determine MSDs above 1 
ns. One must therefore leave open the possibility that the actual 
seff-diffusivity, which by definition is a long-time quantity, may 
. differ from that which is computed over these relatively short 
time scales. Nevertheless, the values reported here appear to 
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Figure 7. Time dependence of In [C^{t)} for [bmim][PF 6 ] at 298 K 
and 0.98 bar. 

be reasonable estimates of the actual self-diffusivity, on the basis 
of the semiquantitative agreement with the results obtained from 
the Stokes -Einstein model. The validity of the self-diffusivity 
can also be tested through use of the following diffusion model. 

Recalling that each ion is surrounded by a cage of other ions, 
we hypothesize that ion diffusion involves a two-step process. 
In the first step, neighboring ions are displaced enough to disrupt 
a cage and form a diffusion pathway. In the second step, an ion 
leaves its cage, moving to another neighboring cage site. To 
test this hypothesis, the characteristic time for the breaking of 
ion cages was determined through use of a cage correlation 
function. 33 The cage for a given ion is determined by the 
neighbor list for that ion, defined as 



l r W>l 



(9) 



where j{nj) is the Heaviside function 



(10) 



and /cm is a cutoff radius for the neighbor list. The cage 
correlation function is defined as 



C^ e (r) = (#(l-«? u '(0,f))>' 



(11) 



where n° u{ is the number of ions that have left ion f s original 
neighbor list at time t and H is the Heaviside function. In 
computing the neighbor lists for all the ions, r cut was set equal 
to the location of the first minimum of the cation/anion radial 
distribution function. Figure 7 shows In C cage (f) for the system 
at 298 K. There is a rapid decay of C cag eM at short times (r < 
0.5 ns), which is attributed to vibrational motion of ions near 
the boundary of a cage. After the period of initial rapid decay, 
QageCO plateaus and then decays exponentially after about 0.8- 
2.0 ns. It is interesting to note that 0.8 ns is about the same 
time in which a shift in the slopes of the MSDs is observed in 
Figure 6. By fitting an exponential function to the region from 
0.8 to 2.0 ns, a time constant of 1.53 ns is obtained for the 
decay of the ion cages. This time constant is indicative of the 
average time needed for any particular ion to leave another ion's 
neighbor list. 

To investigate the rate at which anions exit the ion cages, 
CcageCO was computed by considering only the anions surround- 
ing a given cation. A plot of In C cage (r) at 298 K computed in 
this way is qualitatively similar to Figure 7 and is not shown 



here. The time constant obtained by fitting the intermediate 
region to an exponential is 3.46 ± 0.04 ns. Likewise, the rate 
at which cations exit the ion cages can be studied by computing 
In Ccage(r) by considering only the cations surrounding a given 
anion. The time constant computed for cation departure at 298 
K is- 3.25 ± 0.03 ns. It is not surprising to see that the time 
constant computed from Figure 7 is about half the time constants 
for cation and anion departures. This is because whenever a 
cation departs from an anion's neighbor list that cation 
simultaneously sees the anion depart from its own neighbor list. 
Thus, two cages decay whenever one ion moves out of another 
ion's cage. 

The cage correlation results suggest that diffusion in ionic 
liquids can be modeled semiquantitatively as an activated 
hopping process. Each ion is assumed to reside on a three- 
dimensional lattice, with lattice spacings A given by the distance 
between peaks in g(r). Ions execute a random walk between 
sites with a characteristic time r, assumed to be equal to the 
time constant for cage decomposion. The self-diffusion coef- 
ficient for this model is simply 34 



D 



(12) 



To compute the self -diffusion coefficient for [bmim] using this 
model, we set A = 8 A and r = 3.25 ns. The resulting Aeif for 
[bmim] is 3.28 x 10" 11 m 2 s" 1 . Using A = 8 A and r = 3.49 
ns for [PF 6 ] gives a self-diffusivity of 3.08 x 10" 11 m 2 s" 1 . 
These results are roughly three times higher than the computed 
self-diffusivities and a factor of 10 greater than that estimated 
from the Stokes -Einstein model. It should be expected that this 
simple model over predicts the self-diffusivity, because it does 
not account for correlated hoping motion. In the real system, it 
is likely that an ion that leaves a cage will frequently return to 
its original location, thus decreasing the overall displacement. 
The simplified model does not account for this but rather 
assumes that ions successfully thermalize in a new cage after 
each hop. Nevertheless, this simple model appears to capture 
much of the dynamics responsible for diffusion in this system. 

4.7. Test of Sampling. Figure 6 shows that at 298 K the 
average ion moves only about 2 A per ns of simulation time. 
This sluggish dynamics raises a serious concern over the level 
of phase space sampling in the current simulations, as well as 
the suitability of molecular dynamics for these calculations. As 
a simple test to determine if the simulations were stuck in a 
local potential energy minimum, three additional 300 molecule 
simulations were conducted for [brrumlfPFg] at 298 K. Each 
simulation was started from a different initial configuration and 
at different initial densities.. Figure 8 shows the molar volumes 
of each of these simulations versus time. All three independent 
simulations converged to the same molar volume as that found 
from the original simulation started near the experimental density 
witriin approximately 200 ps. This indicates that, although phase 
space sampling is still a concern for this system, it appears that 
the present calculations have been run long enough to properly 
sample the equilibrium state of the system. 

4.8. Other Dynamic Properties. -The rotational dynamics 
of the anion were investigated by computing the cosine of the 
average angle 0 between the vector from the P atom to the Fi 
atom of an anion, at time zero and time t. The decorrelation of 
this angle with time gives insight into the rotational motion of 
the anion. The rotational time constant was obtained by fitting 
an exponential function to long-time decay of (cos 6). The 
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Figure 8. Plot of the cell volumes of three independent simulations 
of [bmim][PF 6 J at 298 K and 0.98 bar. The horizontal dashed line is 
the average volume of the original simulation. 

computed time constant for anion rotation at 298 K is ap- 
proximately 28.8 ps. In a similar manner, the rotational dynamics 
of the cation were investigated using a vector normal to the 
imidazolium ring. As expected, the rotational time constant is 
much longer at approximately 4.3 ns. This huge separation of 
time scales is significant; although it appears that the transla- 
tional motion of cations and anions is highly correlated and slow, 
their rotational motions occur over vastly different time scales. 

5. Conclusions 

Results of a molecular dynamics simulation of [brnimlfPFe] 
are reported. An all- atom force field for the ionic liquid is 
developed using a combination of ab initio calculations and 
CHARMM22 parameters. The agreement between the experi- 
mental and computed IR spectra is very good, and the vibrational 
motions associated with various' peaks in the experimental 
spectrum are identified. The agreement between experimental 
and computed values of the volume expansivity and isothermal 
compressibility are good, and the agreement between molar 
volumes is excellent. The force field was not adjusted to match 
the experimental data. Discrepancies between the simulation and 
experiment may be due to the use of unoptirnized potential 
parameters or the neglect of polarizability. Liquid structure is 
reported in the form of center-of-mass radial distribution 
functions for cation— cation, cation— anion, and anion —anion 
pairs, as well as site— site RDFs. It is observed that the anion 
tends to orient near the C2 carbon of the cation. Self-diffusion 
coefficients for |>mim] and [PFg] are computed from the slopes 
of the center-of-mass mean-square displacements of the cation 
and anion, respectively. The reported self-diffusivities are 2 
orders of magnitude smaller than the self-diffusivity of water 
at room temperature. The mechanism for diffusion of the ions 
is investigated via the computation of cage correlation functions. 
A simple random walk, diffusion model based on this time 
constant yields a self-diffusivity that is in fair agreement with 
the calculations, as does the Stokes -Einstein estimate based 
on scaling with the viscosity and diffusivity of water. Rotational 
dynamics of the cation and anion are investigated via the 
computation of a common order parameter. The rotational time 
constants are indicative of the very slow rotational dynamics 
of the [bmim] cation but the relatively fast rotational motion of 
the [PFfi] anion. 
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